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ABSTRACT 


Equations of motion characterizing the three dimensional motion of a tethered 
satellite during the retrieval phase are studied. The mathematical model involves 
an arbitrary number of point masses connected by weightless cords. Motion occurs 
in a gravity gradient field. The formulation presented accounts for general functions 
describing support point motion, rate of tether retrieval, and arbitrary forces applied 
to the point masses. The matrix oriented program language MATLAB is used to 
produce an efficient vectorized formulation for a) computing natural frequencies and 
mode shapes for small oscillations about the static equilibrium configuration; b) for 
integrating the nonlinear differential equations governing large amplitude motions. 
An example of time response pertaining to the skip rope effect is investigated. 
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INTRODUCTION AND OBJECTIVES 


Tethered satellites employed in space vehicles have many important practical 
applications. The time response occurring during the retrieval phase is of special 
interest since this motion can sometimes become unstable as the tether length is re- 
duced to zero. One type of motion currently being studied is known as the skiprope 
mode. In this configuration a satellite is initially positioned on the line from the 
orbiter to the earth center (z-axis) and the tether is spinning about this axis in 
a sinusoidal curve resembling the deflection pattern of a skyrope. Some disagree- 
ment currently exists among analysts regarding the dynamics which occurs during 
retrieval. The author’s efforts were directed toward developing a concise computer 
formulation to facilitate a response analysis in this problem. 
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NONLINEAR EQUATIONS FOR LARGE AMPLITUDE MOTIONS 


Consider a tether model described as a series of point masses connected by 
weightless cords assumed to remain in tension throughout the motion. The geometry 
is shown in Figure 1. The vector from mass mj + j to rrij is denoted by r y Thus 

Tj = l j e j (0 j ,(p j ) (1) 

where lj denotes the vector length and (Oj><Pj) are spherical coordinate angles for 
unit vectors ej which specify the member directions. We will assume that lj are 
known functions of time. This allows for treatment of the tether retrieval problem 
during an interval over which l n is reduced to zero at a specified rate. Once a 
particular link has been completely retrieved, the problem dimensionality car be 
reduced by one and the retrieval process can be continued. 

The equations of motion can be integrated numerically when 6 and (p values are 
known for each link. The numerical procedure e below is as follows. 

1. For known 0,<p,6 and <p values, link tensions are computed and these are used 
to compute global accelerations 

2. The global accelerations are used to compute r for each link. 

3. The r values are used to compute 6 and (p values. 

For convenient reference let us call the tether section between m 1+ i and m x as 
link i. The tension force in link i is representable as 

Ti = a, r, (2) 

where a; is a scalar multiplier. A free body diagram of m, is shown in Figure 2. The 
force F{ describes the effects of external loads such as gravity, aerodynamic drag, 
and electromagnetic phenomena. The position of mass i in a rotating reference 
frame tangent to the circular orbit is 

n 

R, = E r i + Rn + X ’ *- n 

j-i 

The length of link i can be computed from 

If = r, • r, = (Ri +1 - Ri ) • (R, +1 - R { ) (4) 
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Figure 1: General Tether Configuration 


Differentiation of this constraint gives 

iiU — (^i+1 “ Ri) * {Rx+l — Ri) (5) 

and 

h “F hh — (^t+i — ^t) * (Ai+i — Ri) + (R \+ 1 ~ Ri) • (Ri + 1 — Ri) (6) 

which is equivalent to 

if + hi i = — ri ■ (iZ t ‘+i - Rt) + r, • r{ (7) 

When a link has constant length then the constraint simplifies to 

r i * (Ri+1 ” Ri) = f t ■ (8) 

We w ri te 

r t = h ei(9i,<Pi) (9) 

so 

■ • de 

r, = l, e, 4- Pi ) = + U(eieO, + e,v Pi) (10) 

I'or the coordinate system of interest e,, and e t # form an orthogonal (but not 
orthonormal triad). Then 
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( 11 ) 


r, • f, = | li e, + U{e x e 8 X + e, v <p,) | 2 


= i? + in\ti*\ 2 % + \e*\ 2 'P 2 i) ( 12 ) 

So 

M. = 1} (|e,5| 2 Of + le.vl 2 ‘Pi) ~ r « ‘ (^«+i ~ #>) ( 13 ) 

- r, • (R i+1 - Ri ) = W* - /-(I «tf| 2 % + M* <P?) = P. ( 14 ) 

The Newtonian equation of motion for a typical mass implies 

R, = — [Fi + cti - 1 r,_i - o,Tj] (15) 

77lj 

so the constraint equation leads to the following relation among link tensions 


r, • (—5— [iq+1 + a, T , - a, + ir i+ i] - — [F t + a,-i r;_i - air,]) = -pi (16) 
m i+ i mi 

— -(r,_i • r,)ai_i + ( — + — )/•«,- — (n • r, + i)a, +1 = p, 

m; m,* m,+i m t +i 

+ r, • [fi/m, - F + i/m, + i] (17) 

Solution of this tridiagonal system determines the member tensions which cau be 
used to compute global accelerations i? t . These global accelerations allow compu- 
tation of 


f, = ~(R t+1 - R,) (18) 

Since r; = /,«,•, we have 

r, = /, e,- + /,(e,-0 <P«) (19) 

r, = /, &i + 2 /,(e,-fl 0; + ei v y>i) 

+ /«(c,-fl <pi + e i8S 8f + e ivv <p\ + 2e iBtfi <pi 8,) (20) 

The last equation can be rearranged as 

r, = /, e, + l{ ei$ 0, + /, e tv) pi + <7, (21) 

Since the coordinate system of interest satisfies the orthogonality condition 
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then 


(22) 


&i8 — — Ci 8 — 0 


implies 

(r, - qi) ■ e, g = 9 t l,\e, e \ 2 

(23) 

Similarly 

0i = (r,-q i )-e ie /(l i \e ig \ 2 ). 

(24) 

implies 

( h - qi) ■ = <fi li\eie\ 2 

(25) 


‘Pi (*"« — ?>) | ). 

(26) 


These relations form the basis for assembling the equations of motion. The 
coordinate reference frame used in the present study assumes a coordinate transfor- 
mation of the form 


r — ( x ->y, z ) — l (sin(0) cos(<^>), sin(0) sin(^), cos(0)) 

which leads to \e^\ = 1 and |e,^| = sin 2 (0). A set of values defines the 

configuration of link i. The values of /, are constant except for i = n, in which 
case, the time dependence of In , in and In are specified by the chosen retrieval rate 
equation. 
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FORCE EFFECTS IN THE GRAVITY GRADIENT FIELD 


We are interested in the dynamics of a tether attached to an orbiter moving about 
the earth with constant angular speed up in a circular orbit of radius pp. A b alance 
of gravitational and centrifugal force requires an orbital speed up = y9P%/Pp where 
pE denotes the radius of the earth. The motion of objects moving relative to the 
satellite are appropriately referred to a frame which rotates with constant angular 
velocity about the earth center. The coordinate frame is positioned such that i is 
pointed tangent to the orbit in the direction of motion, k is directed toward the 
center of the earth, and j is defined as k x i to give a right-handed system. Then the 
angular velocity of the satellite frame is —kup and the position of a generic point 
on the tether is 


r - i x Q + j yp + k z 0 (27) 

The total acceleration of a point referred to the rotating axes is 

a — Rp+uxr + ux(uxr)-\-2uxr + r (28) 

where Rp is a vector from the earth center to the origin of the rotating frame, and 
r and r are velocity and acceleration vectors measured in the rotating frame. The 
total acceleration of point (zcnS/cb^) * s therefore 


a - (i x 0 + j y 0 + k Zp) + i(u% x 0 - 2up Zp ) -f k(ul(pp - *o) + 2u; 0 ip) (29) 

The gravitational force on a particle of mass m is F g = mg p\ R/\R\ 3 where 
R = — i xp — j yp + k (pp — zp). It is helpful to express F g in a power series and 
neglect terms of order p$ 2 . The binomial theorem gives 


F n — mu. 


-i x 0 - j y 0 + k(p 0 - 2 0 )1 [(1 + — 




(30) 


When O(p 0 2 ) is neglected in the laLst equation, it is found that a particle subjected 
to gravity loading and an additional force i F x + j F y + k F z satisfies 


m i i'o + j yp + k zp 


i F x + j F y + k F z + m i(2u 0 z 0 ) 


+ mj(-ulyp) + mk(3up z 0 - 2u> 0 io) (31) 

This e’quation provides useful information regarding feasible two-dimensional states 
of motion in either the orbit plane (i,k) or the plane transverse to the orbit plane 
It is possible to have yp = 0, F y = 0, and all other quantities dependent only 
on the x 0 and zp coordinates. However, if we is assume that F x = 0 and xp - 0, 
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then the gyroscopic term iuoio produces an acceleration in the xq direction. Con- 
sequently, pure motion in the y,z plane is not possible unless zq is negligible. This 
circumstance can occur when a tether, initially in static equilibrium along the z axis, 
is given a small perturbation in the yz plane. Then the motion component for the 
axial direction is of second order magnitude in comparison to the transverse compo- 
nents. The linearized problem for small transverse motions from static equilibrium 
is worth further discussion. 
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TETHER MOTION FOR SMALL PERTURBATION ABOUT THE Z-AXIS 


The problem for motion of a tether experiencing small oscillations about the 
z-direction in the orbital plane were previously investigated by M. Rheinfurth and 
Z. Galaboff [2]. This paper showed that the static tension for extension along the 
z-axis is 


T=?M r u, 0 2 £(a + l-(^) 2 ) (32) 

where Mt is the mass of a tether having length i and a mass M s attached at the 
free end. The parameter a equals 2 M 3 /Mt- For small transverse oscillations, the 
static tension will not change substantially and motion in the z-direction wPl be 
negligible. It is useful to employ dimensionless coordinates according to £ = zjl 
and r = u> 0 t. Following the analysis of Rheinfurth and Galaboff, the linearized 
equations of motion accounting for both in-plane and out-of-plane motion are 


1.5 


d_ 

dZ 




d 2 x 

dr 2 


(33) 


1.5 


d_ ■ 





(34) 


The boundary conditions at £ = 0 are x = y = 0. Force balance conditions on the 
end mass M 3 gives 


dx _ d 2 x 

dy _(Px 
3 V dr 2 

Assuming modal forms X — A(£) sin(u;r) and y = V r (^’) sin(a;r) leads to 


(35) 

(36) 


with corresponding boundary conditions of 


(37) 

(38) 


X(0) = 0, 3.Y'(1) + u 2 X{\) = 0 


(39) 


F(0) = 0, 3T'(1) + (u 2 - 1)7(1) = 0 (40) 
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It shows from the above equations that the modal functions for in-plane and out- 
of-plane motion are identical and a frequency u> T for the in-plane motion also co 're- 
sponds to an out-of-plane frequency 


U y = + 1 (41) 

Reference [2] treated the problem of in-plane motion in considerable detail and gave 
an exact solution for the case where M s = 0, (a = 0). In that case, the exact 
dimensionless frequencies and modal functions are 

u kx = v /6fc(d-.5), k = 1,2,3,... (42) 

XkiO = P2k- 1(0. 0 < £ < 1 (43) 

where Uk x means the k’th in-plane frequency and P 2 k - 1 is the Legendre polynomial 
of degree 2k — 1. Rhein furth and GalabofF showed that the first frequency equals \/3 
and this modal function is a straight line even for the case of a variable mass distri- 
bution. The general Rayleigh quotient formula from [2] expressed in dimensionless 
form is 


“& = [3 f\a+l-e)[X' k (0} 2 /Ly’(1) + 2 f 1 X 2 k (0 d{ 

l Jo J Jo 


(44) 


when M s /Mj > 1 and n > 3, the mode shapes are closely approximated by 


% sin [tv (n - 1)£] (45) 

so the Rayleigh quotient gives 

ukx » - 1) \[l + 3 MJM t (46) 

This formula is accurate within about 4% for n > 3 and M s /Mj > T 

A concise algorithm to compute frequencies and mode shapes for the general 
case can be obtained by employing PC-MATLAB [3] and finite differences as a 
formulation equivalent to that used in [1]. We want to solve 


[( 1 + a - £ 2 )A"(0] = -AX(0 for a = 2 MJM T and 0 < £ < 1 (47) 

The boundary conditions are 


X(Q) = 0 and 2X'(1) = XX(1) 
This corresponds to a modal solution of the form 


(48) 
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( 49 ) 


X(f,t) = X(£) sin (u>o VL5A t ) 

To express this problem in difference form we take 

D = l/N , Zj - j D , 0 < j < N (50) 

The boundary conditions are specified as 

*0 = 0, (X N -. 2 - 4X N -i + ZX N )D~ l = XX N (51) 

and the differential equation becomes 


j 2 -(1 +a)D~ 2 


[X } -i - 2Xj + X i+1 ] + j [X j+l - X 3 ] = AXj, 


2 < j< N - 1. 


(52) 

This gives an eigenvalue problem AX = AX which is solvable using linear algebra 
operators in MATLAB. Figure 2 shows a function tetfrg which determines frequen- 
cies and mode shapes. The figure also lists dimensionless in-plane frequencies for 
several combinations of mass ratio. Computations were made using npts = 50. The 
first four modal vectors for M 3 /Mt = 0.5 are also shown in Figure 3. 
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MATLAB FUNCTION TO COMPUTE TETHER INPLANE FREQUENCIES 


function [ frqs .modevecs ] = tetfrq(npts,msatovmtet) 

* [ f rqs ,modevecs ] = tetf rq( npts ,msatovmtet ) 

■i Frequencies and modal vectors for a tether with an 
% end mass attached in a gravity gradient field. 

S' nrrf* r« o ... 


npts 


msatomtet 


frqs 


number of segments into which the tether 
length is divided for finite difference 
approximation 

sub-satellite mass divided by the 
tether mass 

frequencies for transverse planar motion 
divided by the circular orbit frequency 


■6 Technical reference: Mario H. Rheinfurth and Zachary 
■s J. Galabof f "Modal Analysis of a Nonuniform String 
-s With End Mass and Variable Tension", NASA Technical 
« Paper 2198, August 1983 

% 

% written by Howard Wilson, July 1989 

'O 


a .7 zeros (npts, npts ) ; d = 1/npts; b= ( 2*msatovmtet+l )/d A 2 
a(i,l) - 2 * ( b - 1 ) ; a(l,2) = -b; a(npts,npts-2) = 1/d; 
a ( npts , npts-1 ) = -4/d; a (npts, npts) = 3/d; 
for j = 2:npts-l, a(j,j-l) = j*(j-l)-b; 

a(l,l) = 2* (b- j* j ) ; a(j,j+l) = j*(j+l)-b; 
end 


[ mode vecs, frqs] = eig(a); 

[frqs, id] = sort(sqrt( 1 . 5*diag(frqs) ) ) ; 

modevecs = modevecs( : , id) ; 

modevecs = [ zeros (1, npts ); modevecs]; 


TETHER INPLANE CIRCULAR FREQUENCY DIVIDED BY 
ORBIT ANGULAR SPEED 

Ratio of (Subsatellite Mass)/(Tether Mass) 


freq 

nmbr 

0 

0 . 5000 

1.0000 

1 

1 .7320 

1.7320 

1 .7320 

2 

4 . 2425 

5.5159 

6.7081 

3 

6 . 7066 

10 .1300 

12 . 7360 

4 

9.1551 

14.9072 

18.8904 

5 

11 . 5785 

19.7213 

25.0646 

6 

13.9503 

24 . 5373 

31.2296 

7 

16 . 2319 

29 .3394 

37.3715 

8 

18 .4152 

34.1181 

43 . 4807 

9 

20.5863 

38 . 8659 

49.5493 

10 

22.8428 

43 . 5767 

55.5703 


1.5000 

2.0000 

2.5000 

3 

1.7320 

1.7320 

1.7320 

1 

7.7274 

8.6294 

9.4466 

10 

14.8847 

16.7576 

18.4400 

19 

22.1464 

24.9748 

27.5111 

29, 

29.4180 

33.1951 

36.5799 

39. 

36.6735 

41.3939 

45.6228 

49. 

43.8992 

49.5574 

54.6256 

59. 

51.0851 

57.6751 

63.5773 

68. 

58.2226 

65.7377 

72.4680 

78. 

65.3038 

73.7366 

81.2882 

88. 


Figure 2: Natural Frequency Analysis 
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transverse direction 
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OF POOR QUALITY 


TETHER MODE SHAPES FOR (SATELLITE MASS)/(TETHER MASS) = 1.0 



Figure 3: Tether Modal Shapes 
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COMPUTER IMPLEMENTATION 

The 3D tether dynamics algorithm developed above was implemented in a MAT* 
LAB program having the following modules. 


module 

function 

tetrun 

Create data, call the 
integrator and print results 

tetdat 

Generate data for a tether 
moving in skip rope mode 

tettim 

Compute time when various 
lumped masses have been retrieved 
to the orbiter 

tetlen 

Compute the time dependent 
tether length and 
corresponding derivatives 

tetforce 

Compute gravity and Coriolis 
forces on the tether 

tettop 

Compute position, velocity and 
acceleration of the point to 
which the tether is retrieved 

tetdef 

Form the nonlinear equations 
of motion of the tether 

odeven 

Routine to compute solutions 
at even time increments 

odeTS 

Seventh and eighth order 
Runge Kutta integrator 
based on Fehlberg formulas 


At present, the program has not been checked fully although initial tension values 
computed for the skip rope are reasonable, angular values obtained are large. This is 
apparently caused by the necessity of dividing linear acceleration quantities by ^sin 9 
to get angular accelerations. Figures [4] and [5] show the angular displacements 
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obtained for 6 and The motion rapidly becomes confined to the orbital plane and 
amplitudes of 6 angles become large. The validity of these results seem questionable. 
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Figure 4 
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phi (degree) 
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% Program for 3D Dynamic Simulation of a Tether in a 

% Gravity Gradient Field. 

% by 

% Howard B. Wilson, Summer 1989 

% The chosen data corresponds to a system which was analyzed by 
% David B. Weaver in "Analysis of Transverse Wire Oscillations 
% During The Tethered Satellite System Electromagnetic Mission", 
% November 24,1987 


% global variables for the tether analysis program 

global nlinkO_ lvecO_ lsumO_ cofexp_ gravity_ 
global omegac_ masvecO_ tretrvO_ 
qlobal nlknow lvecnow masvcnow 


% orbit angular speed in rad/sec 
% The orbit period is 1.473 hours 
% gravity constant in ft/sec A 2 
% skip rope spin rate (rad/sec) . This value 
% is 11.31 times the orbit angular speed. 

% number of links used in the model 
% total tether mass in slugs 
% satellite mass 

% satellite initial z coordinate in feet. 

% This equals 20 kilometers 
% Satellite is initially on the z axis 
% Amplitude of initial deflection curve 
% which is half of a sine curve 
% Parameter defining the spatial period 
% length of the initial defelction curve 
% Start initial motion in the x-z plane 
% Initial spin about z axis 

% Initial time derivatives of theta are zero. 

% Time at which frac*(total_tether_ler.gth) 

% has been retrieved 
% The solution is generated until 
% frac* ( total_tether_length) has been 
% been retrieved 

% Generate starting data 

mt=mtet/nlinkO_; masvecO_=[mt+msat;mt*ones(nlinkO_-l , 1 ) ] ; 

[ lvec0_, thvecO , phvecO , thdvecO ,phdvec0 , zO ,x0 ]=. . . 
tetdat ( nl ink0_ , thdl , thdn , phO , phdO , zl , xl , xsin , omeg ) ; 

% Compute times when successive tether masses reach the 
% orbiter 

[ tretrv0_ , 1 sum0_ , cof exp_ ] =tett im ( frac , 1 vec0_ , tf rac ) ; 

yinit=[ thvecO ( : ) ; phvecO ( : ) ; thdvecO ( : ) ;phdvec0( : ) ] ; 

for k=l:nlinkO_, 

nlknow_=nlinkO -k+1; 


omegac_=l . 1847e-3 ; 

gravity_=32 . 2 ; 
omegaskip=. 01134 ; 

nlink0_=5 ; 
mtet=ll . 43 ; 
msat=37 . 7 ; 
zl=65617 ; 

xl=10 ; 

xsin=286 . 315 ; 
omeg=l ; 
ph0=0 ; 

phdO=omegaskip ; 
thdl=0 ; 
thdn=0 ; 
tf rac=7000 ; 

frac=.9; 
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% tmin=tretrvO_(k) ; tmax=. 98*tretrvO_(k+l) ; h=(tmax-tmin)/50; 
tmin=tretrvO_(k) ; tmax=tretrvO_(k+l )-10; h=(traax-tmin)/50; 
tol=l.e-12; trace=l; 

[tout,yout ]=odeven( 'tetdif ' , train, tmax,h,yinit, tol , trace) ; 
ntp=2 ; nt=length ( tout ) ; 

yinit=lagtrp(tretrvO_(k+l) ,tout(nt-ntp+l:nt) ,yout(nt-ntp+l :nt , r ) ) ' 

tira=[ tim; tout ; tretrvO_(k+l ) ] ; yout=[yout;yinit ' ] ; 

ytep= [ ] ; 

for j=4:-l:l, 

yinit( j*nlknow_)=[ ] ; 

yyyy=[yout ( : , ( j-1 ) *nlknow_+l : j*nlknow_) , zeros (nt+1 ,k-l ) ] ; 
ytep= [ yyyy , ytep ] ; 

end 

zout= [ zout ; ytep ] ; 

end 

save tet2 


function [ lng, th, ph, thd,phd , z ,x]=tetdat(n, thdl , thdn,phO ,phdO , zn , . . . 

xO,xsin,omeg) 

% [ lng , th , ph , thd , phd , z , x ] =tetdat ( n , thdl , thdn , pho , phdO , zn , . . . 

% xO,xsin,omeg) 

% Compute data to define an initial planar tether configuration 


% n 

% thdl, thdn 
% 

% 

% phO,phdO 
% 

% 

% zn 
% xsin 


% 

% xO 
% 


Number of links 

Initial values of the time derivative of theta. 
Theta_dot varies linearly from thdl at the first 
link to thdn at the last link. 

Inital values of phi and the time derivative of 
phi. All links have the same inital values for 
these variable. 

The projection of the tether length on the z axis 
Amplitude of a sine wave which is added to a 
linear taper to create the initial radial deflect- 
ion in the x direction. This deflection has the 
form x = xO*(l - z/zn) + xsin*sin(pi*omeg*z/zn) . 
Transverse deflection value at the free end of the 
tether. 


% lng 
% 

% 

% th, ph 
% thd, phd 
% z , x 

% 


Vector of lengths of the individual tether links. 
These lengths are defined by cords on the chosen 
deflection curve. 

Initial theta and phi values. 

Initial time derivatives of theta and phi . 

Initial z and x coordinates viewed normal to the 
plane of the system. 


nn= ( 0 : n ) ' /n ; 
z=zn*nn; 

x=xO*( l-nn)+xsin*sin(pi*omeg*nn) ; 
zdif=z(2:n+l)-z( l:n) ; 
xdif =x( 2 :n+l ) -x( 1 :n ) ; 
ph=phO*ones ( n , 1 ) ; 
phd=phdO*ones ( n , 1 ) ; 
th=-atan2 ( xdif , zdif ) ; 
thd=thdl+( thdn-thdl ) *nn ( 2 : n+1 ) ; 
lng=sqrt ( xdif . *xdif +zdif . *zdif ) ; 
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function [ tvtrv , ltotl , cof ]=tettim( frac, leng, tf ini ) 
% [ tvtrv, ltotl , cof ]=tettim( frac, leng, tf ini ) 


% This function computes the vector of times at which 
% successive links in the tether reach the orbiter during 
% retrieval according to an exponential length reduction rate. 


leng 


The vector of link lengths. The links are 
indexed with the bottom link as the first 
and the base link as the last. 

The fraction of the total length to be 
retrieved by time tfinl. 

The time at which frac*(total tether length) 
has been taken into the orbiter according to 
the exponential rate law. 

The vector of times at which link(n ),..., 
link(l) reach the orbiter. 

The total tether length leng(l )+. . .+leng( 1 ) 
The exponent in the tether equation 
% Tether_length = ltotl*exp(cof *t ) 

n=length( leng ) ; lngsum=f lipx(cumsum( f lipx(leng( : ) ) ) ) ; 
ltotl=sum( leng ) ; cof=log( 1-f rac ) /tfinl ; 
lngsum=[ lngsum; ( 1-f rac ) *ltotl ] ; 
tvtrv=log ( lngsum/1 totl ) /cof ; 


% frac 


tfinl 


fvtrv 

ltotl 

cof 


function [ lvec , lvecd , lvecdd ]=tetlen(t ) 

% [ len,lend, lendd]=tetlen(t) 

% Compute the vector of lengths and their time derivatives. 

% Only the last link changes during the current time interval. 

% global variables for the tether analysis program 

% global nlinkO_ lvecO_ lsumO_ cofexp_ 

% global omegac_ masvecO_ tretrvO_ 

% global nlknow_ lvecnow_ masvcnow_ 

%lvec=lvecO_; lvecd=zeros ( lvec ) ; 1 vecdd=zeros ( 1 vec ) ; 

%break 

% remaining statements are presently inactive, the tether 
% has constant length 

len=lsumO_*exp(cofexp_*t) ; 
lend=cof exp_*len ; 
lendd=cof exp_*lend ; 
len=len-sum( lvecO_( 1 : nlknow_-l ) ) ; 
lvec=len; lvecd=lend; lvecdd=lendd; 
if nlknow_ > 1 

lvec= [ lvecO_( 1 : nlknow_-l ) ; lvec ] ; 
lvecd= [ zeros ( nlknow_-l , 1 ) ; lvecd ] ; 
lvecdd= [ zeros ( nlknow_-l , 1 ) ; lvecdd ] ; 
end 


function f orc=tetf orce ( t , r , v) 


XXXIII-18 


% Force effects due to gavitational attraction and 
% rotation in a gravity gradient field 

% global variables for the tether analysis program 

% global nlinkO_ lvecO_ lsumO_ cofexp_ gravity_ 

% global omegac_ masvecO_ tretrvO_ 

% global nlknow_ lvecnow_ masvcnow_ 

wsq=omegac_*omegac_ ; w2=2*omegac_; 
fx=w2*v( : ,3) . *masvecO_( 1 :nlknow_) ; 
fy=-wsq*r ( : , 2 ) . *masvecO_( 1 :nlknow_) ; 

fz=( ( 3*wsq ) *r( : , 3 ) -w2*v( : , 1 ) ) . *masvecO_(l :nlknow_) ; 
f orc=[ fx , fy , f z ] ; 

% 'in tetforce, (omegac_,masvec_,forc) ' 

% omegac_, mas vec_, fore 
% pause 

%f orc=[ zeros ( nlknow_, 2 ) , ( gravity_*masvecO_( 1 :nlknow_) ) . *ones (nlknow_, 


function [ rtop , rtopd , rtopdd ]=tettop(t) 

% [ rtop, rtopd, rtopdd ]=tettop(t) 

% Coordinates of the tether top and related derivatives. 

% The top is fixed at present. 

rtop = [ 0 ; 0 ; 0 ] ; 

rtopd =[0; 0; 0]; 

rtopdd= [ 0 ; 0 ; 0 ] ; 


function yt=lagtrp(t ,tdat ,ydat) 

% yt=lagtrp(t, tdat,ydat) 

% Lagrange interpolation of a vector function yt which 
% is a function of scalar time. Interpolation is 
% performed such that at tdat(i) the value of yt is 
% ydat(i,:). Thus, yt is a polynomial of degree n-1 
% where n is the number of components in tdat 
[mrow,mcol ]=size(ydat) ; 
yt=zeros( 1 ,mcol ) ; 
for k=l:mrow 

yt=yt+ydat (k , : )*plag(k,tdat,t) ; 
end 


function y=plag( j ,tdat, t) 

% y=plag ( j , tdat , t ) 

% generates a Lagrange interpolating polynomial 

td=tdat ( : ) ' ; td(j)=[]; 

y=prod ( t-td ) . /prod ( tdat ( j ) -td ) ; 


function zp = tetdif(t,z) 

% zp = tetdif(t,z) 

% Dynamical equations for a lumped mass 3D tether model 
% written by Howard B. Wilson, Summer 1989 
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% global variables for the tether analysis program 

% global nlinkO_ lvecO_ lsuraO_ cofexp_ 

% global omegac_ masvecO_ tretrvO_ 

% global nlknow_ lvecnow_ masvcnow_ 

% t, z 
% nlnow_=5; 

n4 = length ( z ) ; n = round(n4/4); w3 = ones(l,3); 
nl = n+1; n2 = nl+n; n3 = n2+n; zn = zeros(n,l); 

% Obtain local angular coordinates and their derivatives 

th ~ z(l:n); ph = z(nl:n2-l); thd = z(n2:n3-l); phd = z(n3:n4); 
ct - cos(th); st = sin(th); cp = cos(ph); sp = sin(ph); 
zp = zeros ( z ) ; zp(l:n2-l) = z(n2:n4); 

% Lengths of tether links and corresponding derivatives defined 
% by an external function 

[ len , lend , lendd ] = tetlen(t); 

% Reciprocals of mass elements needed later 

masinv = 1 ./masvecO_? miv3 = masinv*w3; 

% Compute unit base vectors in spherical coordinates 
% along with spatial derivatives of those vectors 

ur = [st.*cp,st.*sp,ct] ; 
ut = [ ct . *cp , ct . *sp , -st ] ; 

% size(ur) 

% size(zn) 

up = [ -ur ( :,2) , ur ( : , 1) ,zn] ; 
utt = -ur; 

utp = [-ut( : ,2) , ut ( : ,1) ,zn] ; 

U PP = [-ur( : ,1:2) ,zn] ; 
len3 = len*w3; 
thd3 = thd*w3; 
phd3 = phd*w3; 

% Local Cartesian components of position vectors and 
% velocity components 

r = len3.*ur; 

% [size(len3);size(ut) ;size(thd3 ) ; size (up) ;size(phd3 ) ] 
v == len 3 . * ( ut . *thd3+up . *phd3 ) + lend*w3 . *ur ; 

% Compute base point motion defined by an external function 
[RT,RTD,RTDD] = tettop(t); 

RT = RT ( : ) ' ; RTD = RTD( : ) ' ? RTDD = RTDD(:)'; 


% global coordinates and velocity vectors 

R - ones ( n , 1 ) *RT+f 1 ipx ( cumsum( f lipx ( r ) ) ) ; 
V — ones ( n ; 1 ) *RTD+f lipx ( cumsum( f lipx ( v ) ) ) ; 
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% Nodal forces defined by an external function depending on 
% time, position and velocity 

FRC = tetforce(t,R,V) ; 

% Form the tridiagonal equations to solve for member tensions 
if nlknow_ > 1, 
lsq = len.*len; 

b = lsq. * ( [masinv( 2 :n) ; 0 ]+masinv(l :n) ) ; 
c = r( l:n-l, : ) . *r( 2 :n, : ) ? 
c = -masinv( 2 :n) .*sum(c' ) ' ; 
fovm = [ FRC . *miv3 ( 1 : n , : ) ; RTDD ] ; 

e = sum( (r.*( f ovm( 1 : n , : )-fovm( 2 :n+l , : ) ) )' )'<’ 

e = e+sum( (v.*v) ' ) '-lend. *lend-len. *lendd; 

% Solve for tension multipliers 

alp = trisol( [0;c] ,b,c,e) ; 

% Compute member forces 

tenfrc =[[00 0 ] ; ( alp*w3 ) . *r ] ; 

% Compute global accelerations of point masses 

RDD = ( FRC+tenf rc( 1 : n, : ) -tenfrc ( 2 :n+l, : ) ) . *miv3 ( 1 :n , : ) ; 

RDD = [ RDD ; RTDD ] ; 

% Second derivatives of link radii 
rdd = RDD ( 1 : n , : ) -RDD ( 2 : n+1 , : ) ; 
else 

lsq = len.*len? 
b=lsq*masinv( 1 ) ; 
f ovm=FRC . *miv3 ( 1 , : ) -RTDD ; 
e=sum( (r.*fovm)'); 

e=e+sum( (v.*v)' )-lend.*lend-len.*lendd; 
alp=e/b; 

tenf rc=(alp*w3 ) . *r ; 

RDD=(FRC-tenfrc) .*miv3( 1, : ) ; 
rdd=RDD-RTDD; 


end 

% Compute local velocity dependent acceleration terms 

a = ( utt . *thd3+utp . *phd3 ) . *thd3+(upp. *phd3+utp. *thd3 ) . *phd3 ; 

a = a.*len3; 

Compute the second derivatives of the spherical coordinate 
% angles to complete formation of the equations of motion. 

rdda = rdd-a; 

thdd = sum( ( rdda . *ut ) ' ) ' -2*lend. *thd; 


XXXIII-21 


thdd = thdd./len; 

phdd = sum( ( rdda . *up) ')'•/( st . *st)-2*lend. *phd; 
phdd = phdd./len; 
zp(n2:n4) = [thdd; phdd]; 

% End 


function [tout, yout] = odeven(F, to, tfinal, h, yO , tol, trace) 

% 

% INPUT: 

% F - string containing name of user-supplied problem description. 

% Call: yprime = fun(t,y) where F = 'fun'. 

% t Time (scalar). 

% y Solution column-vector. 

% yprime - Returned derivative column-vector; yprime(i) = dy(i)/dt 

% to - Initial value of t. 

% tfinal- Final value of t. 

% yO - Initial value column-vector. 

% h - Time increment between which solution values are output 

% tol - The desired accuracy. (Default: tol = l.e-3). 

% trace - If nonzero, each step is printed. (Default: trace = 0). 

% 

% OUTPUT: 

% tout - Returned integration time points in increments of h 
% (kolumn-vector ) . 

% yout - Returned solution, one solution column-vector per tout-value. 

% 

tout= ( to :h: tfinal) ' ; ntims=length(tout ) ; ncols=length(yO ) ; 
yout=zeros ( ntims , ncols ) ; yout(l, : )=y0( : ) 
for j=l:ntims-l 

[ ttemp , y temp ] =ode7 8 ( F , tout(j), tout(j+l), yout(j,:)', tol, trace); 

[ nr , nc ] =s i ze ( ytemp ) ; 
yout ( j+1 , : )=ytemp(nr, : ) ; 
end 


function [tout, yout] = ode78(F, to, tfinal, yO, tol, trace) 

% ODE78 Integrates a system of ordinary differential equations using 
% 7th order formulas. 

% 

% [tout, yout] = ode78(F, to, tfinal, yO, tol, trace) 

% 


% INPUT: 

% F 

% 

o. 

% 

% 

% to 

% tfinal- 
% yO 
% tol 
% trace - 
% 

% OUTPUT: 


String containing name of user-supplied problem description. 
Call: yprime = fun(t,y) where F = 'fun', 
t - Time (scalar), 

y - Solution column- vector . 

yprime - Returned derivative column-vector; yprime(i) = dy(i)/dt 
Initial value of t. 

Final value of t. 

Initial value column-vector. 

The desired accuracy. (Default: tol = l.e-6). 

If nonzero, each step is printed. (Default: trace =0). 
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% tout - Returned integration time points ( row-vector ) . 

% yout - Returned solution, one solution column- vector per tout-value 

% 

% The result can be displayed by: plot(tout, yout). 

% Daljeet Singh 

% Dept. Of Electrical Engg., The University of Alabama. 

% 11-24-1988. 

% The Fehlberg coefficients: 

alpha = [ 2 ./ 21 . 1/9 1/6 5/12 .5 5/6 1/6 2/3 1/3 101]'; 
beta = [ [ 2/27 00 0 0000000 00] 

[ 1/36 1/12 000000 00000] 

[ 1/24 0 1/8 0000000000]* 

[ 5/12 0 -25/16 25/16 000000 000] 

[ .05 0 0 .25 .2 00000000] 

[ -25/108 0 0 125/108 -65/27 125/54 000000 0 ] 

[ 31/300 000 61/225 -2/9 13/900 000 000] 

[200 -53/6 704/45 -107/9 67/90 300000 ] 

[ -91/108 0 0 23/108 -976/135 311/54 -19/60 17/6 -1/12 0000 

[2383/4100 0 0 -341/164 4496/1025 -301/82 2133/4100 45/82 45/164 18/41 0 0 
[ 3/205 0000 -6/41 -3/205 -3/41 3/41 6/41 000 

[-1777/4100 0 0 -341/164 4496/1025 -289/82 2193/4100 51/82 33/164 12/41 0 
Chi =[00000 34/105 9/35 9/35 9/280 9/280 0 41/840 41/840]'; 
psi = [1 0 0 0 0 0 0 0 0 0 1-1 -1 ]'; 

pow = 1/8; 

if nargin < 6, trace = 0; end 
if nargin < 5, tol = l.e-6; end 

% Initialization 
t = tO; 

% the following step parameters are used in ODE45 
% hmax = (tfinal - t)/5; 

% hmin = (tfinal - t)/20000; 

% h = (tfinal - t)/100; 

% The following parameters were taken because the integrator has 
% higher order than ODE45. This choice is somewhat subjective, 
hmax = (tfinal - t)/2.5; 
hmin = (tfinal - t)/10000; 
h = (tfinal - t)/50; 

y = yo ( : ) ; 

f = y*zeros( 1,13); 
tout = t; 
yout = y . ' ; 

tau = tol * max( norm(y , 'inf'), 1); 

if trace 
% clc, t, h, y 
clc, t, h 

end 

% The main loop 

while (t < tfinal) & (h >= hmin) 

if t + h > tfinal, h = tfinal - t; end 

% Compute the slopes 
f ( : ,1) = f eval ( F , t , y ) ; 
for j = 1: 12 

f(:,j+l) = feval(F, t+alpha( j ) *h, y+h*f *beta( : , j ) ) ; 

end 
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% Truncation error term 
gammal = h*41/840*f *psi ; 

% Estimate the error and the acceptable error 

delta = norm(gammal , 7 inf 7 ) ; 

tau = tol*max(norm(y , 'inf ' ) ,1-0) ; 

% Update the solution only if the error is acceptable 
if delta <- tau 
t = t + h? 
y = y + h*f*chi; 
tout = [tout; t ] ; 
yout = [yout; y . 7 ] ; 

end 

if trace 

home , t , h ; y 
home , t , h 

end 

% Update the step size 
if delta ~= 0.0 

h = min(hmax, 0 . 8*h* ( tau/delta) A pow) ; 

end 
end ; 

if (t < tfinal) 

disp( 'SINGULARITY LIKELY. 7 ) 
t 

end 
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CONCLUSIONS AND RECOMMENDATIONS 


The computer implementation of the 3D version tether dynamics program should 
be studied further to resolve numerical difficulties encountered when one or more of 
the 0-angles are small. An alternative formulation which may be appropriate for this 
case is to develop the equations of motion based on a different set of Euler angles 
than those presently used in the program. This work will be performed during the 
coming months. 
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